Introduction

femR is an R package that solves PDEs on R relying on the Finite Element Method (FEM). Let \(\Omega \subset \mathbb{R}^d\) be the domain of our interest, femR package solves advection-reaction-diffusion problem of the following form:

\[ \begin{cases} -\mu \Delta f + \boldsymbol{b} \cdot \nabla f + c f = u \qquad & in \ \Omega \\ BC \ conditions \ (Dirichlet \ or \ Neumann) \qquad & on \ \partial \Omega, \end{cases} \] where \(\mu \in \mathbb{R}\), \(\boldsymbol{b} \in \mathbb{R}^d\), \(c \in \mathbb{R}\) are the diffusion coefficient, the transport coefficient and the reaction coefficient, respectively. This vignette illustrates how to use femR to solve a partial differential equation. In particular, the following section explains step by step the solution of a simple Poisson problem.

Solving a Poisson problem

Let \(\Omega\) be the unit square \([0,1]^2\). Consider the following problem defined over \(\Omega\):

\[ \begin{cases} -\Delta f = u &in \ \Omega \\ \quad f = 0 &on \ \partial \Omega \end{cases} \]

where \(u(x,y) = 8\pi^2 \sin( 2\pi x)\sin(2\pi y)\) is the forcing term and \(\partial \Omega\) is the boundary of \(\Omega\) where we have prescribed homogeneous Dirichelet boundary condition. The exact solution of the previous problem is \(u_{ex}(x,y) = \sin(2\pi x) \sin(2 \pi y)\). The following figure shows the exact solution of the problem, on the left hand side, and the triangulation of the unit square domain, on the right hand side.

The following window wraps all the steps needed to solve the problem.

## 1. Defining domain
data("unit_square", package="femR")
mesh <- Mesh(unit_square)
## 2. Defining the solution of the PDE
Vh <- FunctionSpace(mesh, fe_order=2) 
f <- f <- Function(Vh)
## 3. Defining the differential operator
L <- -laplace(f) 
## 4. Defining the forcing term and the Dirichlet boundary conditions as standard R functions
# forcing term
u <- function(points){
    return(8.*pi^2* sin(2.*pi*points[,1])*sin(2.*pi*points[,2])) 
}
# Dirichlet BC
dirichletBC <- function(points){
  return(matrix(0,nrow=nrow(points), ncol=1))
}
## 5. Building the PDE object
pde <- Pde(L, u, dirichletBC)
## 6. computing the discrete solution
pde$solve()
## 7. Plots
plot(f)
fig_h <- contour(f)
subplot(fig_exact %>%hide_colorbar(), fig_h, margin = 0.05) %>%layout(annotations = list(
            list(x=0.155, y=1.05, text = "exact solution",    showarrow = F, 
                 xref='paper', yref='paper'),
            list(x=0.875,  y=1.05, text = "discrete solution", showarrow = F, 
                 xref='paper', yref='paper')))

The following sections explain more in detail the previous chunk of code.

Defining the domain

First, we create a proper domain object exploiting the function create_mesh which takes a R named list as parameter. Such named list should contain geometrical information of the domain:

femR contains the unit_square named list that can be used to perform some tests.

data("unit_square", package="femR")
unit_square <- Mesh(unit_square)
plot(unit_square)

Moreover, femR package provides an utility that reads files storing such geometrical information and provided by third-party software such as FreeFem++. Indeed, read_mesh(filename) reads a .mesh file from FreeFem++ and returns a named list that can be passed to create_mesh function. The following chunk, that is not executed, shows how read_mesh works.

filename <- "path/to/domain.mesh" #domain.mesh stored by the FreeFem++ utility savemesh() 
domain <- read_mesh(filename)
mesh <- Mesh(domain)

Defining the solution of the PDE and the differential operator

First, we define initialize a FunctionSpace object passing the mesh tha we previously build and the finite element order. Note that, the package provide first order and second order finite elements. First order finite element is the default parameter. Then, we define a Function, solution of the PDEs, belonging to the function space Vh.

## solution of the PDE
Vh <- FunctionSpace(unit_square, fe_order=2)
f <- Function(Vh)

Then, we define the differential operator in its strong formulation.

L <- -laplace(f) # L <- -div(I*grad(f))
                 # In general:
                 # L <- (-mu)*laplace(f) + dot(b,grad(f)) + c*f

Note that, according to the well-known identity \(div(\nabla f) = \Delta f\), we could have written L<--div(I*grad(f)), where I is the 2-dimensional identity matrix. Indeed, femR package let you define every term of a generic second-order linear differential operator \(Lf= -div(K\nabla f) + \mathbf{b}\cdot \nabla f + a f\) as follows:

Although, the interface let the user define a generic second-order linear differential operator: -div(K*grad(f)) + ..., at this stage the package has been tested only in the isotropic case, i.e K=I, where I is the two-dimensional identy matrix.

Building the PDE object

We define the forcing term of the differential problem and the Dirichlet boundary condition as a simple R functions. Such functions will be provided as parameters to the object which encapsulates the concept of PDEs.

## forcing term
u <- function(points){
    return(8.*pi^2* sin(2.*pi*points[,1])*sin(2.*pi*points[,2])) 
}

## Dirichlet BC
dirichletBC <- function(points){
  return(matrix(0,nrow=nrow(points), ncol=1))
}

Computing the FE solution

Finally, we build a pde object passing the differential operator L, as first parameter, the forcing term u, as second parameter, the Dirichlet boundary condition dirichletBC, as third parameter:

## create pde
pde <- Pde(L, u, dirichletBC)

We can compute the discrete solution of the Poisson problem calling the solve method:

## solve problem
pde$solve()

Since, we are dealing with a simple differential problem and we know the analytic form of the exact solution, we can compute the \(L^2\) norm of the error exploiting on the mass matrix provided by the pde object:

exact_solution <- function(points){
    return( sin(2.*pi*points[,1])*sin(2.*pi*points[,2]))
}
u_ex <- as.matrix(exact_solution(pde$get_dofs_coordinates()))
error.L2 <- sqrt(sum(pde$get_mass()%*%(u_ex-pde$solution())^2))
cat("L2 error = ", error.L2, "\n")
#> L2 error =  0.0004136347

Graphics

The R package plotly is exploited to plot the computed solution.

## plot solution 
plot(f)
# contour 
contour(f)

Note that, the base::plot function and graphics::contourfunction have both been overloaded and return a plotly object that can be modified as one wishes.

plot(f, colorscale="Jet") %>% layout(scene=list(aspectmode="cube"))